Written by E. Montes (emontesh@usf.edu) and Eduardo Klein (eklein@usb.ve) on Auguts 28, 2020.
This code pulls data from NOAA’s ERDDAP servers and creates time series plots of sea surface temperature (SST) and chlorophyll-a concentration (CHL), and maps showing the latest available data for the selected region.
First, let’s load required libraries
library(readr)
library(rerddap)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(flexdashboard)
## Warning: package 'flexdashboard' was built under R version 3.6.2
library(reshape2)
library(leaflet)
library(ggplot2)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'xts'
## The following object is masked from 'package:leaflet':
##
## addLegend
## The following objects are masked from 'package:dplyr':
##
## first, last
library(dygraphs)
library(plotly)
## Warning: package 'plotly' was built under R version 3.6.2
## Registered S3 method overwritten by 'httr':
## method from
## print.cache_info hoardr
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(mapdata)
## Loading required package: maps
library(RColorBrewer)
palette(brewer.pal(8, "Set2"))
Query SST data from ERDDAP
## remove all spaces from string
NoSpaces = function(x){
return(gsub(" ", "", x))
}
## set site coordinates and time for SST extraction
SSTSiteName = "Patagonia" ## for the resulting file name
SSTcoords.lon = -64.
SSTcoords.lat = -41.7
SSTstartDate = "2002-06-01"
## set climatological date start-end
SSTclimStartDate = "2002-06-01"
SSTclimEndDate = "2012-12-31"
## set dataset source
SSTsource = info("jplMURSST41")
##
## Get sst
SST <- griddap(SSTsource,
time=c(SSTstartDate, "last"),
longitude = c(SSTcoords.lon,SSTcoords.lon),
latitude = c(SSTcoords.lat,SSTcoords.lat),
fields = "analysed_sst",
fmt = "csv")
SST = SST[,c(1,4)]
names(SST) = c("time", "SST")
## convert time to a Data object
SST$time = as.Date(ymd_hms(SST$time))
Calculate SST climatology
SST.clim = SST %>% filter(time>=ymd(SSTclimStartDate), time<=SSTclimEndDate) %>%
group_by(yDay = yday(time)) %>%
summarise(SST.mean = mean(SST),
SST.median = median(SST),
SST.sd = sd(SST),
SST.q5 = quantile(SST, 0.05),
SST.q10 = quantile(SST, 0.10),
SST.q25 = quantile(SST, 0.25),
SST.q75 = quantile(SST, 0.75),
SST.q90 = quantile(SST, 0.90),
SST.q95 = quantile(SST, 0.95),
SST.min = min(SST),
SST.max = max(SST))
Plot SST time series
SST.xts = as.xts(SST$SST, SST$time)
dygraph(SST.xts,
ylab = "Sea Surface Temperature (Deg C)") %>%
dySeries("V1", label ="SST (Deg C)", color = "steelblue") %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesBackgroundAlpha = 0.2,
hideOnMouseOut = FALSE) %>%
dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>%
dyRangeSelector(dateWindow = c(max(SST$time) - years(5), max(SST$time)))
## subset SST for last year
SST.lastyear = SST %>% filter(year(time)==max(year(time)))
## make the plot
pp = ggplot(SST.clim, aes(yDay, SST.mean))
pp = pp + geom_line() + geom_smooth(span=0.25, se=FALSE, colour="steelblue") +
geom_ribbon(aes(ymin=SST.q25, ymax=SST.q75), fill="steelblue", alpha=0.5) +
geom_line(data=SST.lastyear, aes(yday(time), SST), colour="red") +
ylab("Sea Surface Temperature (Deg C)") + xlab("Day of the Year") +
theme_bw(base_size = 9)
ggplotly(pp) %>% plotly::config(displayModeBar = F)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Save SST time series data
write_csv(SST, path = paste0(NoSpaces(SSTSiteName), "_SST.csv"))
write_csv(SST.clim, path = paste0(NoSpaces(SSTSiteName), "_Climatology.csv"))
Create a map of the latest SST data
sstInfo <- info('jplMURSST41')
# get latest 3-day composite sst
GHRSST <- griddap(sstInfo, latitude = c(-60., -20.), longitude = c(-90., -47.), time = c('last','last'), fields = 'analysed_sst')
mycolor <- colors$temperature
w <- map_data("worldHires", ylim = c(-60., -20.), xlim = c(-90., -47.))
ggplot(data = GHRSST$data, aes(x = lon, y = lat, fill = analysed_sst)) +
geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
geom_raster(interpolate = FALSE) +
scale_fill_gradientn(colours = mycolor, na.value = NA) +
theme_bw() + ylab("latitude") + xlab("longitude") +
coord_fixed(1.3, xlim = c(-90., -47.), ylim = c(-60., -20.)) + ggtitle("Latest daily SST data")
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Removed 4925876 rows containing missing values (geom_raster).
Query CHL data from ERDDAP
## remove all spaces from string
NoSpaces = function(x){
return(gsub(" ", "", x))
}
## set site coordinates and time for CHL extraction
CHLSiteName = "Golfo Nuevo" ## for the resulting file name
CHLcoords.lon = -63.8587
CHLcoords.lat = -43.1842
CHLstartDate = "2012-01-01"
## set climatological date start-end
CHLclimStartDate = "2012-01-01"
CHLclimEndDate = "2016-12-31"
## set dataset source
CHLsource = info("erdMH1chla8day")
##
## Get CHL
CHL <- griddap(CHLsource,
time=c(CHLstartDate, "last"),
longitude = c(CHLcoords.lon,CHLcoords.lon),
latitude = c(CHLcoords.lat,CHLcoords.lat),
fields = "chlorophyll", fmt = "csv")
CHL = CHL[,c(1,4)]
names(CHL) = c("time", "CHL")
CHL = na.omit(CHL)
## convert time to a Data object
CHL$time = as.Date(ymd_hms(CHL$time))
Calculate CHL climatology
CHL.clim = CHL %>% filter(time>=ymd(CHLclimStartDate), time<=CHLclimEndDate) %>%
group_by(yDay = yday(time)) %>%
summarise(CHL.mean = mean(CHL),
CHL.median = median(CHL),
CHL.sd = sd(CHL),
CHL.q5 = quantile(CHL, 0.05),
CHL.q10 = quantile(CHL, 0.10),
CHL.q25 = quantile(CHL, 0.25),
CHL.q75 = quantile(CHL, 0.75),
CHL.q90 = quantile(CHL, 0.90),
CHL.q95 = quantile(CHL, 0.95),
CHL.min = min(CHL),
CHL.max = max(CHL))
Plot CHL time series
CHL.xts = as.xts(CHL$CHL, CHL$time)
dygraph(CHL.xts,
ylab = "Chlorophyll a (mg m-3)") %>%
dySeries("V1", label ="CHL", color = "steelblue") %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesBackgroundAlpha = 0.2,
hideOnMouseOut = FALSE) %>%
dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>%
dyRangeSelector(dateWindow = c(max(CHL$time) - years(5), max(CHL$time)))
### CHL Last year with smoothed Climatology {data-width=250}
## subset CHL for last year
CHL.lastyear = CHL %>% filter(year(time)==max(year(time)))
## make the plot
pp = ggplot(CHL.clim, aes(yDay, CHL.mean))
pp = pp + geom_line() + geom_smooth(span=0.25, se=FALSE, colour="steelblue") +
geom_ribbon(aes(ymin=CHL.q25, ymax=CHL.q75), fill="steelblue", alpha=0.5) +
geom_line(data=CHL.lastyear, aes(yday(time), CHL), colour="red") +
ylab("Chlorophyll a (mg m-3)") + xlab("Day of the Year") +
theme_bw(base_size = 9)
ggplotly(pp) %>% plotly::config(displayModeBar = F)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#Step 10 Save CHL time series data
write_csv(CHL, path = paste0(NoSpaces(CHLSiteName), "_CHL.csv"))
write_csv(CHL.clim, path = paste0(NoSpaces(CHLSiteName), "_Climatology.csv"))
#Step 11 Create a map of the latest CHL data
require("rerddap")
require("ggplot2")
require("mapdata")
# get latest Monthly chl (VIIRS)
chlaInfo <- info('nesdisVHNSQchlaMonthly')
viirsCHLA <- griddap(chlaInfo, latitude = c(-20., -60.), longitude = c(-90., -47.), time = c('last','last'), fields = 'chlor_a')
# get latest 8-day chl (MODIS)
chlaInfo_8d <- info('erdMH1chla8day')
MODIS_CHLA_8d <- griddap(chlaInfo_8d, latitude = c(-20., -60.), longitude = c(-90., -47.), time = c('last','last'), fields = 'chlorophyll')
# Map monthly chl (VIIRS)
mycolor <- colors$chlorophyll
w <- map_data("worldHires", ylim = c(-60., -20.), xlim = c(-90., -47.))
ggplot(data = viirsCHLA$data, aes(x = lon, y = lat, fill = log(chlor_a))) +
geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
geom_raster(interpolate = FALSE) +
scale_fill_gradientn(colours = mycolor, na.value = NA) +
theme_bw() + ylab("latitude") + xlab("longitude") +
coord_fixed(1.3, xlim = c(-90., -47.), ylim = c(-60., -20.)) + ggtitle("Latest VIIRS Monthly Chla")
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Removed 974803 rows containing missing values (geom_raster).
# Map 8-day chl (MODIS)
mycolor <- colors$chlorophyll
w <- map_data("worldHires", ylim = c(-60., -20.), xlim = c(-90., -47.))
ggplot(data = MODIS_CHLA_8d$data, aes(x = lon, y = lat, fill = log(chlorophyll))) +
geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
geom_raster(interpolate = FALSE) +
scale_fill_gradientn(colours = mycolor, na.value = NA) +
theme_bw() + ylab("latitude") + xlab("longitude") +
coord_fixed(1.3, xlim = c(-90., -47.), ylim = c(-60., -20.)) + ggtitle("Latest MODIS 8-day Chla")
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Removed 677567 rows containing missing values (geom_raster).